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We investigate the simple one-dimensional driven model, the totally asymmetric exclusion process, 
coupled to mutually interactive Langmuir kinetics. This model is motivated by recent studies on 
clustering of motor proteins on microtubules. In the proposed model, the attachment and detach¬ 
ment rates of a particle are modified depending upon the occupancy of neighbouring sites. We first 
obtain continuum mean-field equations and in certain limiting cases obtain analytic solutions. We 
show how mutual interactions increase (decrease) the effects of boundaries on the phase behavior of 
the model. We perform Monte Carlo simulations and demonstrate that our analytical approxima¬ 
tions are in good agreement with the numerics over a wide range of model parameters. We present 
phase diagrams over a selective range of parameters. 

PACS numbers: 87.16.Uv, 05.70.Ln, 87.10.Hk 


I. INTRODUCTION 

Driven diffusive systems show a very rich behavior. 
Even one-dimensional systems exhibit boundary induced 
phase transitions[IfjBj with a complex phase behavior. 
One such model is that of a totally asymmetric exclusion 
process (TASEP)fTfHU] coupled to the Langmuir kinet¬ 
ics (LK)[B]. In that model a single species of particles 
performs unidirectional hopping on a ID lattice. The 
particles are assumed to have hard-core repulsion which 
prevents more than one particle from occupying the same 
lattice site. Such a system is coupled to Langmuir kinet¬ 
ics by allowing for adsorption (desorption) of particles at 
an empty (filled) lattice site with fixed respective kinetic 
rates. It was shown in Refs. 0 0 that the combina¬ 
tion of TASEP and LK results in nonconserved dynamics 
with unusual features such as the appearnce of a high- 
low coexistence phase separated by stable discontinuities 
in the density profile. The novel phase behavior has its 
origin in the competing kinetics of TASEP and LK. How¬ 
ever, in the thermodynamic limit, it is expected that the 
bulk effects are predominant with boundaries becoming 
insignificant. In fact, the competition between bulk and 
boundary dynamics can occur only if one rescales the 
attachment (detachment) kinetic rates 0 Ej such that 
they decrease with increasing system size in a particular 
fashion. 

Besides being fundamentally interesting, understand¬ 
ing nonequilibrium physics of driven systems is of par¬ 
ticular interest in biological systems [Bj |TTj. One such 
particular system is that of molecular motors performing 
directed motion along one-dimensional molecular tracks. 
Typically kinetic rates are such that the fraction of 
track over which the motor moves before detaching is 
finite [T2j. This allows for the bulk dynamics to com¬ 
pete with the boundary, potentially giving rise to unusual 
nonequilibrium stationary states. Recently, exclusion 
process on networks have been used to model cytoskeletal 
transport [TB . It was shown that active transport pro¬ 
cesses spontaneously develop density heterogeneities at 
various scales. An important aspect that needs to be in¬ 


cluded in a study of motor transport is that of mutual 
interactions (MI) between motors. Seitz et al|T3] ob¬ 
served that in presence of an obstacle, a molecular motor 
walking on a microtubule tends to stay attached for a 
longer time. Muto et al[T5] reported on long-range co¬ 
operative binding of kinesin to a microtubule. The de¬ 
tachment could depend on the biochemical state of the 
motor PEI [T7J which might itself be determined by the 
presence or absence of neighbouring motors. 

In a recent study on kinesin-1 motors moving on 
microtubule|18]. the authors performed numerical simu¬ 
lations of binding/unbinding dynamics incorporating mu¬ 
tually attractive interaction between the motor proteins. 
Their results were in agreement with the experimental 
observation, in particular clustering of motors on micro¬ 
tubules. Mututal interactions in addition to the hard¬ 
core repulsion introduce additional correlations as in the 
Katz-Lebowitz-Spohn (KLS) model, which is a generic 
model of interacting driven diffusive systems|T9]. By 
modifying the hopping rate of particles depending upon 
the occupancy of next nearest neighbour, the model gives 
rise to exotic features such as localized downwards shocks 
and phase separation into three distinct regimes[20], 
However, in the case of molecular motors, due to the mu¬ 
tual interactions, the attachment and detachment rates 
of a motor molecule are modified depending upon the 
state of the neighbouring sitesjT8j. Assuming that the 
hopping rate is unaltered, this corresponds to the ordi¬ 
nary TASEP (with no correlations besides the hard-core 
repulsion) with density (local) dependent LK. 

In this paper, we focus on the TASEP coupled to mu¬ 
tually interactive Langmuir kinetics. We investigate how 
mutual interactions can tilt the balance in favor of pre¬ 
dominantly bulk effects by enhancing LK. We show that 
that this is indeed the case when both the attachment 
and detachment rates are enhanced significantly due to 
the mutual interactions. In the case of the kinetic rates 
being significantly reduced due to the interactions, one 
suppresses the bulk effects giving rise to rich and com¬ 
plex phase behavior. We also explore the more interest¬ 
ing scenario in which the kinetic rates are modified in an 
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Figure 1: A graphic representation of the TASEP model 
with Langmuir dynamics. Particles are injected at the first 
site with rate a and extracted at the last site with rate /3 .The 
particles move exclusively to the right, with unitary rate, if 
the next site is empty. In this model the only interaction 
between the particles is the hard-core repulsion, this means 
that hopping over particles and multiple occupation are not 
allowed. The injection, extraction and hopping to the next 
site constitute the TASEP model. The Langmuir dynamics 
consists of the detachment and attachment from and to the 
background. The attachment and detachment rates are oj a 
and lod respectively. For analysis of this model see for ex¬ 
ample Refs. 03 E]. Mutual interactions are incorporated by 
modifying the attachment and detachment rates; see Fig. [2] 


asymmetric fashion. The paper is organized in the fol¬ 
lowing way. In Sec. [TTJ we present the model composed 
of TASEP coupled to the modified Langmuir Kinetics. 
We first present the modfication to the LK due to the 
mutual interactions. We then obtain continuum mean 
field equation describing the steady state density profile 
of particles on the lattice. In Sec. |III| we study three dif¬ 
ferent cases of modified LK. We first consider the case in 
which the unmodified LK rates are assumed to be equal 
and mutual interaction enhances (suppresses) the attach¬ 
ment and detachment rates by the same factor. The sec¬ 
ond case corresponds to the unmodified LK rates being 
equal, but the mutual interaction enhances one and sup¬ 
presses the other by the same factor. The last case is the 
most general one in which all the model parameters are 
freely chosen. We do not explore this case in detail. In 
Sec. |IV| we summarize our findings. 


II. THE MODEL 

The model consists of a ID lattice with sites i = 
1,2, ...,1V; see Fig. [T] Each site on the lattice can ei¬ 
ther be occupied by one particle or no particle. There 
are three different sub-processes that govern the dynam¬ 
ics of the system: 

(a) Particles are injected at site 1 with rate a and ex¬ 
tracted at site N with rate /?. 

(b) Particles at site i = 1, ...N — 1 can hop to site i + 1 
if site i + 1 is unoccupied. 

(c) Particles at site i = 2,..., N — 1 can detach from the 
lattice with rate wjj and attach to site i = 2,..., N— 1 
with rate ui A . 
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Figure 2: The mutual interaction are incorporated in the 
TASEP model with LK dynamics by modifying the attach¬ 
ment and detachment rates if neighbouring sites are occupied. 
For each occupied neighbouring site the attachment rate is 
multiplied by 5 and the detachment rate is multiplied by 7 . 


All the rates are defined such that the hopping rate is uni¬ 
tary. Processes a and b are the dynamics of the TASEP 
model, process c is the interaction with the background. 
The interaction with the background is called Langmuir 
kinetics. We note that the only interaction between the 
particles is assumed to be the hard-core repulsion. 

For the sake of completion we show the equations for 
site-occupancy below. These equations are the same as 
reported in Ref. |5j. The equation for the occupancy of 
each site is given by 

4 n »(*) = rii-i(t)[l - rij(i)] - n,(t)[l - n i+1 (t)\ 

+ w A [l-n i (t)\-u> D rii(t), (1) 

with rii the occupancy at site i, which can either be one 
or zero. The equations for the boundary sites are 

^m(t) = a[l - m(t)] - ni(i)[l - n 2 (t)\, (2) 

^n N {t) = njv-i(t)[l - n N (t)] - /3n N (t). (3) 

The mutual interactions of the particles are included 
in the equation by modifying the attachment and detach¬ 
ment rates to a and 10d respectively. The attachment (de¬ 
tachment) rate if both neighbouring sites are unoccupied 
is uj a (wp). If either the left or right neighbouring site is 
occupied the attachment (detachment) rate becomes 6 uia 
('yajo), and if both neighbouring sites are occupied S 2 uj a 
(7 2 wd); see Fig. 0 In our model, the hopping rate of par¬ 
ticles on the lattice is unaltered in presence of mutual in¬ 
teractions. This is in contrast with the KLS model where 
the hopping rates are modified according to the occu¬ 
pancy of nearest and next-nearest neighbours [19] whereas 
the binding/unbinding kinetics remain unaltered. 

In order to include the mutual interactions of the par¬ 
ticles the following substitutions are needed: 

oja —> kM[l + (ui + i-|-?Tj_i)(<5—l) + 7ii + i7ii_i(<5 2 — 25+1)], 

(4) 
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uj d -» w_ D [l + (n i+ i+n i _i)( 7 -l)+n i+ in i _i( 7 2 - 27 +l)]. 

(5) 

At present it is not clear how already bound mo¬ 
tors modify the binding kinetics of motors. It has 
been suggested that presence of a bound motor could 
change the lattice locally somehow leading to a modified 
binding/unbinding kinetics [2Tj. However, the underlying 
mechanism remains unkown. 

In order to obtain useful solutions for the distribution 
of particles on the lattice, the following two steps are 
required. First one goes from the equation for the occu¬ 
pation of the sites, where each site can have either value 
one or zero, to an equation of the average occupation of 
the sites. Second in the limit of large system sizes a semi- 
continuous variable x instead of the discrete parameter 
i is used for the position on the lattice. This method is 
the same as used in [B], The average density at a site is 
defined as ( m(t )) = Pi(t). In stationary state the average 
(rii) can either be a time or a sample average. In order 
to take the averages of Eqs. 000 with substitutions 
0- ([5] the higher order correlations are needed. Instead 
of solving these equations exactly a mean field approach 
is used which consists of the approximation: 

(rii(t)n i+1 (t)) « {ni(t))(n i+ i(t)). (6) 

The lattice constant e is defined as e = L/N. For 
simplicity the length of the lattice is fixed to one L = 
1. For large system sizes, TV 1 the quasicontinuous 
position variable x = i/N is introduced. This means that 
the average density at site i is now defined as (rij(i)} = 
p(x,t). The equation for the average density profile in 
stationary state, to leading order in e becomes 

0 = - d^p + (2 p — 1 )d x p + ffyjl + p(5 - 1)] 2 (1 — p) 

-H D [l + p( 7 -l)]V (7) 

The [1 + p{8 — l )] 2 and [1 + p(j — l )] 2 parts of the equa¬ 
tion are due to the mutual interactions. In equation [7] 
the total detachment/attachment rates are used, defined 
as flA = ojaN and = uijj TV. The equations for the 

boundary sites (Eqs. 00) become the boundary condi¬ 
tions. 

p(0) = a, p(l) = 1-/3. ( 8 ) 

One can now take the continuous limit, e —> 0 for a nor¬ 
malized lattice this means TV —> 00 . In order to ensure 
that the attachment/detachment rates per unit length do 
not become infinite, the total rates Qa and are kept 
constant. In the continuous limit, e —> 0, the second order 
differential equation |7]) becomes a first order differential 
equation, but the two boundary conditions remain. This 
means that the problem is overdetermined. However one 
can find solutions to the equation in the continuum limit 
that satisfy one of the two boundary conditions. The 
full density profile is is constructed from the possible so¬ 
lutions. The crossover position from one solution to an 


other is obtained by matching the currents j(x) jBj. 

j(x) = p(x)[l - p(x)\ (9) 

For a normalized lattice in the continuous limit the 
crossover region is localized and a discontinuity in the 
density profile appears. Though the crossover region in 
this case is localized it does span a finite number of sites 
implying that in the case of a finite sized lattice the 
crossover region spans a finite fraction of the normalized 
lattice. 

The model without MI exhibits a particle-hole symme¬ 
try, in the sense that a particle attaching to the lattice 
means that a vacancy detaches from the lattice and vice 
versa. The same holds for a particle entering (leaving) 
the system on the first (last) site, which can be seen as 
a vacancy leaving (entering) the system. And a parti¬ 
cle hooping to a neighbouring site on the right equals a 
vacancy hopping to the left. Due to this symmetry the 
following transformations 

rii(t) O- 1 - n N -i(t ), (10a) 

a -H- /?, (10b) 

WA O (10c) 

leave Eqs. 0. © and ([ 3 ]) invariant [6|. However this 
particle-hole symmetry is no longer apparent if MI is in¬ 
cluded. 


III. MODIFIED LANGMUIR KINETICS 

In this section the solutions of Eqn. |7]) in the contin¬ 
uum limit are presented for three different cases. The 
first case corresponds to where LK rates are both en¬ 
hanced or reduced simultaneously by the same amount, 
and second case to where ffy is enhanced while is 
reduced by the same amount. The third case is the most 
general one in which the attachment and detachment 
rates are independently modified. This case ix explored 
in least details. With these solutions the density profiles 
are constructed and compared with Monte Carlo simu¬ 
lations of the model. Besides the density profiles phase 
diagrams are made which show the characteristics of the 
solutions for different values of a and /?. 

Case 1: mutual interaction with enhanced LK rates 

The model simplifies significantly in the case that 
£Ia = Qd = Q and 8 = 7 = 1 + r). This means that 
the attachment and detachment rates are both multiplied 
by 1 + 17 for each occupied neighbouring site. Values of 
77 < — 1 result in negative rates, therefore 77 is restricted 
to values larger than -1. Positive 77 increases and neg¬ 
ative r] decreases the LK dynamics if neighbouring sites 
are occupied. In this case Eqn. ([7| in the continuous limit 
becomes 


0 = (2 p- 1 ) ( d x p - H [1 + prj\ 2 ) . 


( 11 ) 
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The special case where 77 = 0, the case without mutual 
interactions, corresponds to the symmetric case analysed 
in [Sj. Eqn. (11 1 has three solutions. A constant solution 
pi = 1/2 which is the Langmuir isotherm. The other 
solutions are p a and pp which obey the left and right 
boundary condition respectively. 


a + (1 + r/a)flx 
1 — (1 + rja)r]£lx 


( 12 ) 


_ i - P + Mi - P) + iMs -1) n o\ 
1 - ( 77(1 - (3) + 1)t)Q,{x - 1) 

The full density profile p(x) is a combination of these 
solutions 


Pa 

for 0 < X < x a , 


Pi 

for x a < x < xp, 

(14) 

pp 

for xp < x < 1. 



Where x a and xp are obtained by equating the currents 
of the solutions, j a (x a ) = ji and jp(xp) = ji |S|. The 
domain of the solutions can be explained as an attraction 
of the density to the Langmuir isotherm as one moves 
away from the boundary. As one moves away from the 
boundary the influence of the bulk dynamics i.e. the 
Langmuir kinetics becomes more dominant and therefore 
the density tries to reach the Langmuir isotherm pi. 

In the case that x a > xp the constant solution pi is 
not part of the density profile and the profile becomes 


phase consists of the p a solution on the left side and the 
pp on the right side of the lattice. The phase diagrams 
are constructed using the information of the domain of 
each of the solutions. The detailed construction of the 
phase diagram with 77 = 0 is reported in Ref. jB]. 

The phase diagrams for nonzero 77 are shown in Fig. [3] 
The behaviour of the phase diagrams for changes in 77 
are similar to changing O. For increasing 12 the area 
occupied by the LD, LD-HD and the HD phase in the 
phase diagram decrease and and eventually disappear . 
The key difference between changing Q and 77 is that the 
influence of 77 is stronger in the phases containing the 
HD phase. As seen in Fig. [3]T) and (d), for increasing 
77 the HD phase disappears quickly from the phase dia¬ 
gram, while the LD phases decreases slowly. This is due 
to the high probability of occupied neighbouring sites in 
regions of high density. Changes in 12 and 77 do not ef¬ 
fect the area of the MC phase, which is always confined 
to the upper right quarter of the phase diagram. If one 
keeps increasing 77 eventually only the LD-MC, MC, MC- 
HD and LD-MC-HD phase remain in the phase diagram. 
Further increasing 77 does not change the phase diagram 
any more, however the density profiles do change. For 
77 —» 00 the LD and HD parts of the density profile oc¬ 
cupy an infinitely small domain on the boundaries of the 
profile. This means that due to the increase in Langmuir 
dynamics the density on the whole lattice is equal to the 
Langmuir density pi. 


Density profiles 


p a for0<x<x w , ^ 

pp for x w < x < 1, 

where x w is obtained by matching the currents of the two 
solutions, j a {x w ) = jp{x w ) [Bj. In this case a discontinu¬ 
ity appears at x w . 

Phase diagrams 

There are seven characteristic density profiles, called 
phases. Depending on 12 and 77 all or some can be present 
in the phase diagram. We follow the same terminology 
as in Ref. (BJ. The simplest three are the high density 
(HD), the low density (LD) and the maximum current 
(MC) phase. In the HD (LD) phase the density is higher 
(lower) than 1 / 2 , and the density profile is given by the 
Pf 3 {pa) solution. In the MC phase the density profile is 
equal to 1/2 over the whole domain. This is called the 
maximum current phase because the current is maximal 
for p = 1/2. The density profile in the MC phase does 
not depend on the boundary conditions a and f3. 

Due to the interaction with the background it is pos¬ 
sible for two or all three solutions to coexist in a density 
profile. This happens in the LD-HD, LD-MC, MC-HD 
and the LD-MC-HD phases. For example the LD-HD 


With equations for p a , pp and pi the density profiles 
are constructed. In Fig. [4] the density profiles for the five 
different phases with ST = 0.3 and 77 = 2.0 are shown 
,these correspond to the phase diagram in Fig. [3](d). The 
first thing to notice is that the p a and pp solutions are not 
straight lines, in contrast to the solutions for 77 = 0 which 
are straight lines. For 77 > 0 the p a and pp solutions are 
concave up with a positive slope. This can be explained 
by an increase in attraction to the Langmuir isotherm 
as the density increases. For example in figure [4] (a), 
if one moves away from the left boundary the density 
increases. This increase in density leads to an increase in 
mutual interactions. In the case of positive 77 this results 
in an increase in LK dynamics and therefore an increase 
in attraction to the Langmuir isotherm. This increase 
in attraction causes the slope to increase. The p a and 
pp solutions with 77 < 0 are concave down with positive 
slope, this can be explained with the same arguments as 
in the case of 77 > 0 . 

The analytical solutions of the density profiles in the 
continuum limit are compared with Monte Carlo simula¬ 
tions of the model; see Fig. [4] Due the the finite size of 
the lattice used in the simulations one can expect certain 
discrepancies between the analytical results and the sim¬ 
ulations. If the p a or the pp solutions are not part of the 
density profile, one or both of the boundary conditions 
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Figure 3: Five different phase diagrams for S2 = 0.3 and 
different values of 77 : (a) r/ = 0, (b) rj = 0.5, (c) 77 = —0.5, 
(d) r/ = 2.0 and (e) 77 = —0.9. (a) Corresponds to the case 
without mutual interaction. The influence of Q is analysed in 
[ 6 j. From (b) and (d) it becomes clear that due to the increase 
in 77 the HD region disappears more quickly from the phase 
diagram than the the LD phase. This is explained by the fact 
that the mutual interactions are most apparent in region of 
high density. I 11 cases (c) and (e) 77 is decreased. Again it 
becomes clear from these figures that a change in 77 has more 
influence on the phases containing HD regions than on phases 
containing LD regions. 


are not met. This happens in the LD, LD-MC, MC and 
the MC-HD phase. In these phases a so called boundary 
layer forms where the boundary condition is not met j.6|. 
A boundary layer is a discrepancy between the analytical 
result of the equation in the continuum limit (Eqn. @) 
and the simulation results of the model with a finite sized 
lattice. This discrepancy occurs at the boundary where 
the boundary condition is not met, and will occupy an 
increasingly small domain for an increasing number of 
lattice sites used in the simulations. One can also expect 
a discrepancy between the analytical density profile and 
the simulation results where there is a crossover from one 
solution to the other. In the analytical density profile the 
crossover is localized on the scale of the rescaled variable 
x. However if the lattice has a finite number of sites, 
the crossover will span a finite fraction of the normalized 
lattice. 


Figure 4: Density profiles for = 0.3 and 77 = 2, this 
corresponds to the phase diagram in Fig. [ 3 ] (d). The boundary 
conditions are (a) a = 0.01 /3 = 0.01, (b) a = 0.2 (3 = 0.2, (c) 
a = 0.8 f3 = 0 . 1 , (d) a = 0.8 /3 = 0.8 and (e) a = 0.1/3 = 0.8. 
The solid lines are the analytical solutions for the density 
profiles (Eqs. PPl Q, the constant solution pi is included 
as a dash-dot line. The solid lines with noise are the result 
of the Monte Carlo simulations with a lattice of 1000 sites, 
averaged over 2000 simulations. The analytical solutions for 
fl = 0.3 and the same boundary conditions, but without MI 
(77 = 0) are included as a dotted line to emphasize the effect 
of the mutual interactions, (a) The MI do change the density 
profile significantly, but do not change the phase, (b) Due to 
the MI the LI< dynamics increases, and the profile changes 
from a LD-HD phase for 77 = 0 to a LD-MC-HD phase, (c) 
The MI cause a phase change from the HD phase to the MC- 
HD phase, due to the increased LK dynamics.(d) There is no 
difference between the solution with or without MI. (e) The 
MI induce a phase change from the LD phase to the LD-MC 
phase, due to the increased LK dynamics. There are two 
types of discrepancies between the analytical results and the 
Monte Carlo simulations. Boundary layers are formed at the 
left boundary in (c) at the right boundary in (e) and at both 
boundaries in (d). Other discrepancies between the analytical 
result and the simulations occur where there is a transition 
between the p a , the pp and or the pi solution. Both types of 
discrepancies can be explained by the finite size of the lattice 
used in the simulations. 


Case 2: mutual interactions with antisymmetric 
modified LK rates 

In the previous case the MI increased the LK dynamics, 
which is not the attractive interaction as reported in ra¬ 
in this case the model is analyzed for attractive interac¬ 
tions. Again Qa and fare set equal, = O. 

But in this case the mutual attraction are incorporated 
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in an antisymmetric manner, S is increased and 7 is de¬ 
creased by a factor ip, 5 = 1 + ip and 7 = 1 — ip, with 
— 1 < ip < 1. Depending on whether ip is positive or 
negative the interactions between the particles is respec¬ 
tively attractive or repulsive. In this case Eqn. Q in the 
continuous limit becomes 


0 = (2p-l)d x p+fil[l +pip] 2 (l-p)pp’] 2 p. (16) 

In order to find solutions for the density profile, the equa¬ 
tion is simplified by neglecting the terms of order ip 2 . In 
the next section the limit of this approximation is dis¬ 
cussed. With this approximation Eqn. (16 1 simplifies to 

0 = (2p-l)d x p+n~2n(l-ip)p. (17) 


This equation has the same form as the one for the model 
without MI but with unequal f l a and Qd (Eqn. ® 0 , 

0 = (2 p — 1 )d x p + — (Qd + £Ia)p- (18) 


Eqn. ( |18[ ) exhibits a particle-hole symmetry 
(Eqn. (10)j7 therefore this symmetry is also appar¬ 
ent in Eqn. (Jl7|. However this is not a property of the 
model and is only apparent if the ip 2 terms in Eqn. (16) 
are neglected. The transformation 


p{x) ->• 1 - p{l - x), 

(19a) 

a 

(19b) 

P —■ ► a, 

(19c) 

n ->• 0(1 - 2ip), 

(19d) 

1 , 

(19e) 

leaves Eqn. (fl7| invariant. Using this transformation 
density profiles for —1 < ip < 0 can be obtained from so- 


lutions to Eqn. ( fl7| with 0 < ip. Therefore the analysis is 
restricted to positive values of ip. Solutions to Eqn. ( fl 8 j ), 
obtained by [ 6 ], are Lambert W functions. Using these 
solutions one finds that the solutions to Eqn. ( fl7| for 
ip > 0 are 


= 2<r7o {w -' H,wl +11 + s 


for a < 1 / 2 , 

( 20 ) 


PP 


2 (T=t 0 (Wo [y(x)] + 1 ) + \ for 1 - /? > p h 
2 ( 1 -^) (Wo [~y(x)} + 1) + \ for \ < 1 - /3 < pi. 


( 2 !) 

where p a obeys the left and pp the right boundary con¬ 
dition. y(x) Is given by 


y{x) = 

2Q 


1 — ip 


exp 


ip 

(1-VO 2 

ip 


{2po - 1 ) - 1 


(x - x 0 ) + 1 ^ ( 2 p 0 - 1 ) - 1 


( 22 ) 


with po = a, Xg = 0 for p a and po = 1 — p3, Xq = 1 for 
pp. The constant solution pi = * s the equivalent 


of the Langmuir isotherm in the case without MI. The 
density profile is "attracted" to this constant solution, 
as explained in the previous section. The solution p a is 
stable only for a < 1 / 2 , and pp for pi < 1/2 [B]. 

The full density profile is constructed from the solu¬ 
tions obeying the left and right boundary conditions, and 
calculating x w , the position of the transition from p a to 
pp, by matching the currents of these solutions. 


Phase diagrams 

Using the same method as in the previous case the 
phase diagrams are constructed from the information 
about the domain of the solutions. There are four pos¬ 
sible phases, these have the same characteristics as the 
phases of the model without MI but with 7 ^ ftp [B], 
Due to the similarity only a short explanation is given 
here, for a more elaborate discussion of the phases one 
can consult |5-. In the LD (HD) phase the full density 
profile is governed by p a (pp); boundary layers appear 
at the right (left) end of the lattice. The condition for 
the LD phase is x w > 1 and a < 1/2, and for the HD 
phase the condition is x w < 0 and pi < 1/2. The M phase 
occurs for /3 > 1/2 and x w < 0. This phase is called the 
"Meissner" phase due to similarities with the Meissner 
phase in super conducting materials | 6 ] . In the M phase 
the density profile is independent of both boundary con¬ 
ditions, therefore boundary layers occur at both ends of 
the lattice. Because the solution is not stable for these 
values of pi the profile is given by pp{\/2) [ 6 j|. This phase 
can be seen as the equivalent of the MC phase in case 
without MI |U or case 1, because a profile in the MC 
phase is also independent of the boundary conditions. 
The LD-HD phase, where phase coexistence occurs, is 
split in two regions. In the region (3 < 1/2 the profile 
obeys both the left and right boundary condition. In 
the region p3 > 1/2 only the left boundary condition is 
obeyed. The right part of the density profile is inde¬ 
pendent of the right boundary condition and is given by 
pp{l/2), this phase is indicated as LD-HD(M). Profiles in 
the LD-HD (M) phase have a boundary layer at the right 
end of the lattice. The conditions for the LD-HD phase 
are 0 < x w < 1, a < 1/2 and f3 < 1/2. For the LD- 
HD (M) phase the conditions are 0<x IU <l, a < 1/2 
and f3 > 1/2. In Fig. [ 5 ] three phase diagrams are shown 
for different values of ip. 

Because Eqs. ( |20| and © are derived using the ap¬ 
proximation ip 2 = 0, the phase diagrams are also an ap¬ 
proximation which hold in the limit of small ip. 


Density profiles 


Using Eqs. (20) and (21 1 the density profiles can be 
constructed. The domain of each of the solutions is 
determined by matching the currents of the solutions, 
Eqn. ([£]). In Fig. [ 6 ] five density profiles are depicted, one 
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Figure 5: Three phase diagrams obtained with Eqs. (201 
and @. With Q. = 0.3 and (a) ip = 0.001, (b) ip = 0.4, (c) 
ip = 0.8. The phases which contain the p a solution disappear 
from the phase diagram for increasing MI, until only the M 
and HD phases are left. The area of the M phase in the phase 
diagram occupies an increasingly large part of the upper half 
of the phase diagram for increasing ip. This stands in contrast 
with the MC phase in case 1, where the MC phase is confined 
to the upper right quarter of the phase diagram and the area 
is independent of any parameter; see Fig. [3] 


each for a phase in phase diagram [5] (b) (fl = 0.3, ip = 
0.4). It is clear from the Fig. [6] that there is good agree¬ 
ment between the simulations and the analytical solu¬ 
tions (Eqs. ( 20|21 )). The apparent discrepancies between 
the analytical and numerical results are caused by the fi¬ 
nite size of the lattice used in the simulations. Besides 
this there is also a small discrepancy between the an¬ 
alytical solution and the simulations caused by the ap¬ 
proximation ip 2 = 0. This can cause a discrepancy in the 
domain wall position, as can be seen in Fig. and (d). 


Approximation limits 


In deriving Eqs. (201 and (211 terms of the order 


ip 2 were neglected. The neglected part of Eqn. (161 is 


f hp 2 p 2 — 2ttip 2 p 3 , which shows that every ip 2 is coupled 
to either a p 2 or a p 3 . This means that the break down 
of the approximation is governed by both p and ip and 
the approximation hold for low densities regardless of the 
value of 0, because MI do not play a significant role in 
low densities due to the low probability of having occu¬ 
pied neighbouring sites. Fig. [7] illustrates some of the 
limits of the approximation. From Figs. |6](b),(d) and [ 7 ] 
(b) it becomes clear that the domain of the low and high 
density solution can differ significantly even if the p a and 
pp differ only little from the simulation result. 



(b) 
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Figure 6: Density profiles for fl = 0.3 and ip = 0.4, this 
corresponds to the phase diagram in Fig. [5] b. The injec¬ 
tion/extraction rates are (a) a =0.01, /? = 0.7 (LD phase), 
(b) a = 0.1, /3 = 0.7 (LD-HD(M) phase), (c) a = 0.8, /3 = 0.8 
(M phase), (d) a = 0.05, ft = 0.1 (LD-HD phase) and (e) 
a = 0.3, P = 0.2 (HD phase). The dashed lines are the 
Lambert W function solutions to Eqn. (17 ) obeying the left 
and right boundary conditions. The parts of these solutions 
that make up the density profile are represented as solid black 
lines. The dash-dot line is the constant solution. The ana¬ 
lytical solutions for the same values but without MI are in¬ 
cluded as dotted lines to emphasize the influence of the MI on 
the density profiles. Solid lines with noise are the results of 
the Monte simulations with a lattice of 1000 sites, averaged 
over 2000 simulations. Over all there is good agreement be¬ 
tween the simulations and the analytical result. There are two 
causes for the discrepancies, the finite size of the lattice and 
the approximation ip 2 « 0. Due to the finite size of the lat¬ 
tice boundary layers are formed at the right end of (a),(b),(d) 
and at the left end of (c) and (d); and the domain walls in 
(b) and (d) are not localized. The discrepancies caused by 
the approximation are visible in two ways, the density profile 
does not fully coincide with the analytical result and due to 
this the position of the domain wall is shifted. This is visible 
in (b) and (d). 


Case 3: mutual interactions with arbitrarily 
modified LK rates 


Until now, we have considered enhancement or sup¬ 
pression of LK rates in a symmetric or antisymmetric 
fashion. The most general case in which all the relevant 
parameters (f Ia, Qd, <5, 7 , a , /3) are assigned randomly 
chosen values is extremely difficult to treat analytically. 
Due to the large parameter-space ( 6 -dimensional), one 
cannot gain much insight by performing numerical simu- 











































Figure 7: The same legend as in Fig. [6] is used. In order 
to illustrate the limits of the approximation ip 2 ss 0 used in 
deriving the equations for the density profiles, four extreme 
cases are shown. For all figures D = 0.3 was used an for (a) 
ip = -0.99, a = 0.1, 0 = 0.1, (b) ip = 0.8, a = 0.1, 0 = 0.9, 
(c) ip = 0.5, a = 0.9, 0 = 0.1 and for (d) ip = 0.9, a = 
0.9, (3 = 0.1. From (a) it is clear that the approximation holds 
for low densities, regardless of the values of ip. The analytical 
solution in (b) does not fully coincide with the simulations. 
There are two causes for the discrepancies. First of all there 
is a boundary layer on the right side caused by the finite size 
of the lattice. Secondly the pp solution is higher than the re¬ 
sults from the simulation, this is caused by neglecting the ip 2 
terms. Over all the simulations are in good agreement with 
the analytical solutions, however the domain of the solutions 
differs significantly. The analytical profile is in the HD phase, 
x w < 0. The simulation result, on the other hand, is in the 
LD-HD phase, 0 < x w < 1. (c) For this value of ip there 
is agreement between the analytical solution and the simula¬ 
tions, even though density is high, (d) The analytical result 
does not coincide with the simulations due to the combination 
of high density and a ip close to 1. In this case the density of 
the analytical result exceeds 1, which is physically impossible. 


lations. Therefore, we have not explored the most general 
case in any details. However, we consider few represen¬ 
tative cases in which the choice of model parameters is 
based on the observation that the resulting density pro¬ 
files have interesting features when contrasted with the 
case with no mutual interactions. In Fig. [8] we show 
profiles corresponding to four different sets of parame¬ 
ters. As can be seen the profiles look very different from 
those with no mutual interaction. However, as mentioned 
above, at present our analysis of the most general case is 
very qualitative and highly superficial. It is obvious that 
it requires much further investigation. We leave detailed 
analysis of the general case to a future study. 




X X 

Figure 8 : The density profiles for (a) Oa = 0.5, Q.d = 1, 
S = 2, 7 = 0.1, a = 1 and 0=1, (b) Qa = 0.5, Dd = 0.5, 
8 = 2, 7 = 0.5, a = 0 and 0 = 1, (c) Da = 0.5, Q. D = 0.5, 
5 = 2, 7 = 0.5, a = 0 and 0 = 0, (d) D A = 0.5, D d = 1, 
8 = 3, 7 = 0.1, a = 0 and 0 = 0. Solid line with noise 
is the result from the simulation, the dashed line is the so¬ 
lution without MI. In all cases the attachment/detachment 
is increased/decreased, which results in a higher density. In 
(b) and (c) the density overlaps with a small part of the ana¬ 
lytical solution without MI. This can be explained by the in 
low density regions the effect of the MI is small. Though in 
(d) the solution without MI the density is low, the lattice is 
almost completely filled due to the increase/decrease in at¬ 
tachment/detachment. All parameters in (a) and (d) are the 
same except for the boundary conditions, which prevents the 
lattice in (a) from filling up completely. 


IV. CONCLUSION 

In this work, we investigate the simple one-dimensional 
driven model- the totally asymmetric exclusion process, 
coupled to a modified form of Langmuir kinetics. This 
model is motivated by recent studies on clustering of mo¬ 
tor proteins on microtubules. Without addressing the un¬ 
derlying mechanism, it is assumed that the attachment 
and detachment rates of a particle depend on the occu¬ 
pancy of the nearest neighbours of any given site. Ig¬ 
noring density correlations, we obtain continuum mean- 
field equation describing the density profile on the lat¬ 
tice. Imposing certain conditions, we obtain analytical 
solution to the equation and demonstrate using Monte 
Carlo simulations that our analytical solutions are accu¬ 
rate over a wide range of parameters. We show that when 
both attachment and detachment rates are enhanced due 
to mutual interactions, bulk-effects start dominating the 
phase behavior of the model. The two-phase coexistence 
(low and high density) observed in absence of mutual in¬ 
teractions can become three-phase coexistence (low and 
high density with maximum current phase) when mutual 
interactions are attractive. On varying the mutual in¬ 
teraction between particles (attractive or repulsive), we 
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obtain a very rich phase diagram describing the behav- Langmuir kinetics, 
ior of driven diffusive system with mutually interactive 
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